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C ^ Abstract. We show that existing Runge-Kutta methods for ordinary differential 

^J , equations (odes) can be modified to solve stochastic differential equations (sdes) 

with strong solutions provided that appropriate changes are made to the way 
stcpsizcs are selected. The order of the resulting sde scheme is half the order 
of the ode scheme. Specifically, we show that an explicit 9th order Runge-Kutta 
method (with an embedded 8th order method) for odes yields an order 4.5 method 
for sdes which can be implemented with variable stepsizes. This method is tested 
" ,..' ■ by solving systems of sdes originating from stochastic wave equations arising from 
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1. Introduction 

Stochastic wave equations play an important role in the quantum theory of 
decoherence and measurement P E] as well as in computational many-body physics 
0I|S1|S]. Solutions of master equations for completely positive dynamical semigroups [5] 
and for Redfield theory[2] can be expressed as expectations of diadics formed from 
wavefunctions obeying stochastic wave equations. Recently it has been shown that 
exact solutions of the N-body Schrodinger or Liouville-von Neumann equations can 
be expressed as averages of Hartree products of single-body wavefunctions or densities 
which obey stochastic wave equations PIE El El- Such methods could have important 
applications in chemistry and condensed matter physics. Stochastic differential 
equations (sdes) are also widely employed in other areas of physics, engineering and 
finance 13 . 

Unfortunately, efficient numerical techniques for solving such equations have not 
yet been developed. Algorithms in the literature have not substantially improved on 
the primitive methods described ten years ago in the well known text by Kloeden and 
Platen[S]. Methods applicable to general systems of stochastic differential equations 
with multiple Wiener processes have not exceeded an order of 2. The low order of 
such methods restricts their domain of usefulness to one or few equations, and the 
larger systems of equations of interest in physics cannot be solved. 

Recently, one of us noted||9] that with minor modifications classical methods 
for ordinary differential equations (odes) can be used to solve sdes with strong 
solutions. The technique was demonstrated by solving a wide range of low dimensional 
sdes with known exact solutionsj^J- Here we expand upon this idea by developing 
variable stepsize (i.e. adaptive) explicit Runge-Kutta based integrators for sdes. We 
demonstrate the use of the method by solving a variety of stochastic wave equations 
arising in decoherence problems TO, and in stochastic decomposition of the many-body 
problemOJlElEl. 

2. Stochastic Taylor expansion 

There is a close connection between Taylor expansions of solutions of odes and Taylor 
expansions of strong solutions of sdesip. Consider a finite set of sdes, 

m 

dXl = a^ {Xt , i) dt + y IP^ (Xf , t) dW^ , (1) 



fe=i 



represented in Ito 8 form, where j = 1, . . . , n. Here Xt = {X} , . . . , X") and the dWi 
are independent and normally distributed stochastic differentials with zero mean and 
variance dt (i.e. sampled A^(0, dt)). The stochastic variables Wi are Wiener processes. 
Assume that the coefficients a^ and &^ have regularity properties which guarantee 
strong solutions, i.e. that Xl are some fixed functions of the Wiener processes, and 
that they are differentiable to high order. We may then view the solutions of QJ as 
functions Xl = Xj{t, W^, . . . , W™) of time and the Wiener processes. The solutions 
can therefore be expanded in Taylor series. Keeping terms of order dt or less then 
gives 

X' ^xi + ^dt+y^ dw^ 

k=l ^ 
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fc,Z— 1 ^ ^ 

The product of differentials dW^dW} is equivalent to S^^idt in the It6,8, formulation 
of stochastic calculus, so that 



dX't+dt - ^hdt - ^/ - [^ + 2 E i^] ^^ 



k=l 



Comparison to QJ allows us to identify the first derivatives 

dxl _ ^, 



y^dw,^ (3) 



a# = ^^^(^-*) (4) 



^t 



^ m 



,2 v:* 



d^X: 



fc=l i=l * 



(5) 



From these first order derivatives, expressed in terms of a^ and &^, higher order 
derivatives can be computed. Thus a Taylor expansion of the solutions 

X't+At = XI + ^At + y 1^ AW,^ 

t + At t Q^ /_^ Qy^k t 

k—1 '■ 

L y /;^« , Aw,'Awl + ... (6) 



2 ,^ dWfdW. 

k,l=l ^ '■ 

can be obtained for finite displacements At and AW^. 

3. Runge-Kutta methods for sdes 

This Taylor expansion of strong solutions of sdes can be employed to develop Runge- 
Kutta algorithms and other integration schemes[5]. As an example consider the classic 
fourth order Runge-Kutta scheme with four stages 

Kj -/,(X,,,f,) 



2 ' ' 2 



Kf -/,(X,, + -K2,t, + -Ai) 
Kf -/,(X,, +K3,i,+i) 

Xt.^, -Xi, + i(Ki + 2K2 + 2K3 + K4) (7) 

6 



with fj (Xt , t) defined via 



f (^Xt,t) = ^At + y ^AW, 

k=l '■ 



k 
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^t 



+ J2bUX,,t)AWt. (8) 

fc=i 

Here ti is the initial time and ti+i — ti + At. Taylor expansion shows that X^^^^ differs 
from the exact solution by terms of order higher than At^ (i.e. terms of higher order 
than At^, At{AWt^)^, (A^^/=)^ (AW^^)^ (AWl)^ , and (AWf )2AW^/AW7). Thus, this 
stochastic Runge-Kutta algorithm is very similar to its classical counterpart except 
that its order is 2 not 4. 

Generalizations to higher order Runge-Kutta schemes are straightforward. One 
simply replaces the usual stage evaluations of Runge-Kutta with evaluations of 0. 
Since is order At^/^ rather than order At, the order of the sde method is half 
that of the ode method. SDE methods of order 2 and 4 constructed in this fashion 
have been shown to be very accurate in fixed stepsize calculations for small systems 
of sdesPJ- In this manuscript we adapt a 9th order Runge-Kutta methodfT^ with 16 
stages into an order 4.5 method for sdes. 

However, fixed stepsize Runge-Kutta methods are neither accurate nor efficient 
for general systems of equations. To solve the large systems of sdes that arise in 
physical problems we need some means of controlling the local error. 

4. Adaptive stepsizes 

Local error is typically controlled in Runge-Kutta schemes for odes via the use of 
embedded lower order methods [TUl [TTl [T^ [T^ . That is, Runge-Kutta methods can 
often be found wherein a method of order I with k > I stages has an embedded Runge- 
Kutta scheme of order I — 1 which uses some subset of the k stages of the higher 
order method. Differences in the two solutions can be compared to a user requested 
tolerance to decide whether a contemplated step can be accepted or whether a smaller 
stepsize should be considered. Thus local error can be estimated, and stepsizes adapted 
to ensure the accuracy of the solution, at negligible extra cost. Well implemented 
examples of this approach are the 5(4) and 8(7) embedded pairs of Dormand and 
Prince (see ^21 and JS|) respectively) which form the basis of the ode software package 
RKSUITE. A 9(8) pair has been derived by Tsitouras^l] although this algorithm has 
not been included in any software of which we are aware. Runge-Kutta methods of 
order 10 are known^3] but embedded lower order pairs have not been reported. 

Variable stepsize one step schemes such as Runge-Kutta are popular because 
they are simple to understand and easy to implement. Multi-step schemes, such 
as Predictor-Corrector^, which store and use information from previous steps are 
however often much faster and more accurate. Unfortunately, it is not clear how the 
stochastic Taylor expansion developed above can be incorporated into a multi-step 
scheme. Predictor- Corrector J^, for example, employs Lagrange type interpolation 
formulae (to fit fj (Xj , i) at a set of times) , which are explicitly integrated over a 
time interval, to construct both the predictor and the corrector. It is far from clear 
how an analogous scheme would work for the m + 1 variables t, Wi, . . . , Wm- Thus, 
Runge-Kutta methods seem to be the easiest to develop for sdes. 

Once an error has been judged too large to be acceptable, ode codes simply try 
a smaller step and all information about the original step is lost. This procedure 
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obviously cannot yield unbaised solutions in an analogous scheme for sdes. Therefore 
measures must be taken to ensure that the original Wiener process is maintained. One 
way of doing this is to halve the original step and to generate stochastic differentials 
on the two subintervals such that their sum is the original step. This approach was 
originally proposed by Gaines and Lvons[T7| and is known as the method of binary 
Brownian trees. More sophisticated strategies have since been developed 18^ but they 
are specific to individual algorithms and cannot be easily adapted for our purposes. 
A number of schemes have been proposed for choosing the stochastic differentials on 
the subintervals |17l [TO| . The correct approach appears to be that of LambafT^ who 
generates the stochastic differential on the first subinterval by sampling the conditional 
probability 

- -^-oo '^^a.dXbp{Xa)p{xb)5{lS.Wa - Xa)5{Xa + Xfc - AVF) 

/^ dXadXbp{Xa)p{xb)S{Xa + Xb ~- AW) 

- ' exp{-i^^%^4^} (9) 



^27r(At/4) 2(Ai/4) 

where AW is the original stochastic differential, AWa is the stochastic differential on 

the first subinterval and 

1 x^ 

p(x)^^^ cxp{- „,^ ,A (10) 

^2{At/2) ^ 2(At/2)^ ^ ' 

is the independent density for differentials on the subintervals. This implies that the 
stochastic differential AWa on the first subinterval has mean AW/2 and variance 
At/ A. The stochastic differential on the second subinterval AWb must then be given 
by AWb = AW — AWa in order to maintain the original Wiener process. 

Thus, Runge-Kutta methods can be developed for sdes with strong solutions from 
Runge-Kutta methods for odes, and a binary tree variable stepsize strategy can be 
implemented with sampling on the subintervals via Lamba's method 19 . To show 
that the combined approach yields an accurate numerical method we solve a variety 
of stochastic wave equations from the recent physics literature. With the adaptive 
stepsize strategy we have chosen there is a good correlation between the speed of an 
algorithm and its order. We thus chose to employ the highest order Runge-Kutta 
pair available, which as far as we are aware is the 9(8) pair of Tsitouras_14_. We 
have sucessfuUy used other lower order methods such as those of |12l and ^^, but 
for consistency all results reported in this paper were calculated using the method 
described in jTlj . 

5. Examples 

Here we solve 3 sets of equations from the recent physics literature using the order 4.5 
Runge-Kutta method implemented with variable stepsizes as described above. 

5.1. 

The first example we consider is the Gisin-Percival^ stochastic wave equation for the 
nonlinear absorber (Eq. 4.2 of Ref. ^) 

d\^{t,Wt)) = .l{a^ - a)\^it,Wt))dt 

+ {2{a^^)a^ - a^^a^ ~ \{a^)\^)\^{t,Wt))dt 

+ V2{a^ - {a^))\^{t,Wt))dWt (11) 
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where a denotes the usual harmonic oscillator lowering operator. We chose an initial 
state |\l/(0,0)) = |0) where |0) is the lowest eigenstate of a'^a. We choose Wq = in 
this and later examples and the Wiener process Wt is real. The notation (F) indicates 
the quantum expectation (^|y|^). The ensemble average over statistical realisations 
of the Wiener processes is denoted via M\Y] for any Y . The quantity of interest for 
this example is the average density 

p{t)=Km{t,Wt)){^>{t,Wt)\] (12) 

which obeys the deterministic master equation 

^^ = .Ifflt _ a, pit)] + 2a2p(t)at2 _ a^^a^pit) - p{t)a^^a\ (13) 
at 




300 



Figure 1. Mean occupation number nt vs. time t 



To implement our approach we need to find the derivatives of |^(i, T4^()) with 
respect to t and Wt. We immediately see that 



S^^.V2ia'-W))m.«.)) 



(14) 



and using © we also determine that 

^IMII) = .i(,t _ ,)|^(,, Wt)) + ((a^) - {aY)\^{t^ Wt)) 

~ (a' - {a'))'\^it,Wt)) + ((atV) - \{a')\')\^{t, Wt)) 

+ {2{a^')a' ^ a^'a' ~ \{a')\')\^it,Wt)). (15) 

From these results we can now construct fj (Xt , t) using Eq. (jHJ . The Runge-Kutta 
scheme can thus be implemented as discussed above. The dynamics was solved in 
a basis consisting of the lowest 11 eigenstates \n) of a''a with n = 0, . . . , 10. Thus, 
including real and imaginary parts of (n|^(i, Wt)) for n = 0, . . . , 10 our equations 
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consist of a total of 22 real nonlinear coupled stochastic equations. A relative 
tolerance of 10"^'^ was requested. In Fig. 1 we plot the mean occupation number 
Ut = Ti{a''ap{t)} vs time for 50000 stochastic realisations (dashed curve) and for an 
exact solution of Eq. (|13|l performed in the same basis set (solid curve). Agreement 
is very good. 

Due to the increasing number of Wiener processes in the following examples, we 
drop the Wt's as arguments of the stochastic wavef unctions. This change of notation 
is necessary but regretable since this dependence is an essential requirement of the 
method we are testing. 

5.2. 

The second example is the Gisin-Percival^ stochastic wave equation for a quantum 
cascade with absorption and stimulated emission (Eq. 4.4 of Ref. 0) 

d|*(t)> = - .li{a'' + a)\^{t))dt 

+ (2(ata) a^a - {a^af - {{a^a)f)\-^{t))dt 

+ .01(2(a^)a - a^a - \{a)\^)\^{t))dt 

+ V2{a^a-{a^a))\^(t))dWt^ 

+ .lV2{a - {a))\^>{t))dW^. (16) 

Here again we chose the initial state |^(0)) = |0). There are now 2 real statistically 
independent Wiener processes (i.e. M[dWldWf] = 0). The quantity of interest is 
again the density (|12|l which in this case obeys the master equation 

-^ = _ .li[at + a, p{t)] + 2atap(t)a^a - (a^a)V(i) ~ p{t){a''af 

+ .02ap{t)a^ - .Olatap(i) - .01p{t)a^a. (17) 

The derivatives of |^(t)) with respect to t, W^ and Wf are given by 

4^=V2(ata-(ata))|*(t)) (18) 



dW^ 



^Ml ^.lV2(a-(a))|vI.(i)) (19) 

«^^-.H(aUa)l*(t)) 
+ {2{a^a) a^a - {a^af - {{a^a)f)\^{t)) 
+ .01(2(at)a-ata-|(a)|2)|*(i)) 

- (flta _ (ata))2|*(i)) + 2{{{a^a)^) - {a^ a)^)\^ (t)) 

- m{a - {a)r\^{t)) + .01((ata)) - \{a)\^)\^{t)) 

+ .Ol{{a')-{ar)\^{t)). (20) 

The same basis set and tolerance as in example 1 were employed. 

Again we calculated the mean occupation number for 50000 trajectories (dashed 
curve) and for an exact solution of Eq. 1)171) (solid curve). These quantities are 
plotted in Fig. |2Jl. Agreement is again good but convergence is somewhat slower 
than in example 1 since we now have twice as many Wiener processes. 
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Figure 2. Mean occupation number nt vs. time t 



5.3. 

The third example consists of stochastic wave equations for a stochastic decomposition 
of the Schrodinger equation for HehumjH]. 

Neglecting nuclear motion about the center of mass, the Helium wavefunction 
$(ri,r2,i) obeys the deterministic Schrodinger equation (in atomic units h — 1, 
rUe = 1, and e = 1) 



d<^{ri,r2,t) 
dt 



i7^2$(ri,r2,t) 



1. 



V^ 



2 

ri 



2 

r2 



1 



(21) 
-}$(ri,r2,t), 



2 ' 2 - ri r2 |ri-r2| 

for any specified anti-symmetric initial state <i>(ri,r2, 0). Here ri and r2 denote 
positions of electrons 1 and 2 with respect to the nucleus. For our example calculation 
we choose an initial wavefunction of the form 

$(ri,r2,0)=/3(*i(ri,0)*2(r2,0)-*2(ri,0)*i(r2,0)) 

(where (3 = Xj \J2{\ — |(\l'i(0)|^2(0))p) is a normalization factor) which is obviously 
antisymmetric in ri and v-2,- Note that we are implicitly incorporating the two- 
component electron spins into the definitions of ^i and ^2- For our purposes it 
is important that (^i(0)|^2(0)) ^ 0. The actual initial conditions for this example 
calculation were chosen randomly as a mixture of Is and 2s He+ states for each 
electron. 

It can be shown "B" that the exact deterministic wavefunction $(ri,r2,t) evolving 
from (|5.3I) can be decomposed into stochastic waves via an average of the form 

$(ri,r2,t)=/3Af[^i(ri,t)^2(r2,i)-*2(ri,t)*i(r2,t)] (22) 
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where *I/i and ^2 satisfy stochastic wave equations 

12 ^ 



2 

p 

2 



+ ^^a>,(0,)i(0,)2*i(r,t)]di 



s=l 



+ ^V^iZJKO, -(0,)i)*i(r,i)dM^/ (23) 



s=l 



It' ^' 2Re{(vI.,|*2)} *^(^'^)'^ 

12 ^ 

d*2(r,i) = H(--V2--)*2(r,i)-z^c^.(Os)iOs*2(r,t) 

. p 
^^c.,(0,)2(0,)i*2(r,t)]di 



^ V^zZJ(0, - (0,)2)*2(r,t)dW^/ (24) 



P /iTf^iiTi_\r//nt/n \_ \in \_121 



^' ^' 2Re{(vI'i|*2)} '^ ^ 

We have used a notation (F)j ~ {'^ j\F\^> j) in the above equations. Here the tOg and 
operators Os arise through the one-body expansion of the Coulomb interaction 

1 ^ 

—=^c.,0. (1)0,(2) (25) 

In -1-21 ^^^ 

which we performed numerically in a basis of He"^ eigenstates 6' . In the calculation 
reported here p = 8 which means that there are eight real Wiener processes. Since 
the initial states are of s type we included only the basis functions of s type with a 
principle He^ quantum number of 4 or less[n|. This means that the total number of 
equations was 32. Clearly this is by far the most computationally difficult of the three 
examples. 

The derivatives of ^j(r, t) are given by 






d^j{r,t) 
dt 



= V^i^{Os-{0,),)^j{r,t) 
t 

p 



= -*(-Jv2 - -)^,{v,t) -iY^u^s{Os)kOs^,{v,t) 



s=l 



-^c.,(0,)i(0,)2VI',(r,t) 

s=l 



^' ^' 2Re{(vI/i|vI/2)} '^^ ' 



2.=i 
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\J2WMolo,), ~\{o,),\'')^,{v,t) 



10 
(26) 



where k ^ j and j,k — 1,2. From these equations we can now construct fjCKiji) 
using Eq. JHJ). 




Figure 3. Re ('I'{0)|'I'{t)) vs. t 




Figure 4. Im (<I'(0)|1'(t)) vs. t 
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A numerical problem arises in Eqs. H25|) when the overlap (^i|^2) becomes 
small. The terms inversely proportional to this factor vary rapidly and the speed of 
integration slows down greatly. This occurs every few atomic units. Fortunately, this 
can be easily avoided by adding a small piece of ^2(1", t) to ^i(r, i), renormalizing the 
wavefunctions, and carrying the new norm as a weight factor in the stochastic average. 
The antisymmetric nature of the full wavefunction guarantees that this manipulation 
makes no change in the solution. 

In Figs. © and |J3J we plot the real and imaginary parts of (<i>(0)|$(t)) for 200000 
trajectories (dashed curve) and for the exact solution (solid curve) of the Shrodinger 
equation (|21() . Agreement is satisfactory with some deterioration of accuracy as time 
proceeds. 




Figure 5. He energy spectrum 



Finally, we computed the energy spectrum via 



I{E)^ 



nh 



Re 



(vl/(0)|vl/(t))exp 



iEt 



dt~ (*(0)|(5(-B-H2)|*(0))(27) 



which we plot in Fig. ((S)). Again satisfactory agreement is obtained. 

The fact that complete convergence is not achieved even with 200000 realisations 
may be due to the relatively large number of Wiener processes. Unfortunately, 
examples with very large numbers of Wiener processes will arise when the method 
of stochastic wave equations described in [S| is applied to larger atoms or molecules. 
Thus it may be necessary to explore some form of importance sampling to improve 
convergence for these simulation methods. 

6. Discussion 



The numerical strategy discussed in this manuscript provides a method for solving the 
large sets of coupled nonlinear sdes which arise in physical problems. This is currently 
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the best strategy for solving systems of sdes like those that arise from stochastic wave 
equations. However, the large number of stages (16 for JJ) required for high order 
Runge-Kutta formulae limit the efficiency of our approach. Multistep methods for odes 
such as Predictor-Corrector 16 typically require only two evaluations of derivatives per 
step and can be implemented to any desired order. If such methods could be adapted 
for sdes the gain in efficiency could be enormous. Unfortunately, to implement such 
a strategy would require interpolation in m + \ variables with variable stepsizes (here 
m is the number of Wiener processes). At present we do not see how this can be 
accomplished. 
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